% T1IRExecute: Estimates parameters and errorbars for T1 recovery in an IR sequence 
% model = A - B*exp(-TI/T1)

% Two ROIs need to be drawn: First a ROI with only noise and second the ROI
% where the mapping has to be performed

% Output are dicom images of the R1map, R1errormap, T1map and T1errormap, saved
% in the directory of the first image of the serie. IMPORTANT: The values of the R2 and
% R2 errormap are multiplied with 10000 to get into sensible gray values.
% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.

% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007.

clear all

roibased = 0; % 0: pixel based mapping. 1: roibased mapping.
niftioutput = 0; % 0: only .dcm output. 1: .dcm and .nii output


studydir='C:\Users\Henk Smit\Desktop\MRData\110329 Jasper dGemric\dGEMRIC_data\dGEMRIC_patient_001\dGEMRIC1\dgemric001_3420181_1858\T1';
file='0024';

%maskdir='C:\Users\hsmit\Desktop\MRData\Jasper\dGEMRIC_data\dGEMRIC_patient_002\Registered\Mappings_dG1';
%maskfile='0025_T1map';
%maskinfo = dicominfo([maskdir,'\',maskfile]);
%maskimage = uint16(dicomread(maskinfo));

outputdir='C:\Users\Henk Smit\Desktop\MRData\110926 Jasper 10 dGe Volunteers\qmap_1_test_1_ufuk\Mappings';

[file,chosendir]=uigetfile('*.dcm');
cd(chosendir);
cd ..
studydir=pwd;

cd(studydir)
d=dir;
nfol=length(d);
fitdata(1).folder=char(d(1+2).name);

%get the dicom info
for i=1:nfol-2   
    fitdata(i).folder=char(d(i+2).name); 
    info = dicominfo([studydir,'\',fitdata(i).folder,'\',file]);
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime;
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
    fitdata(i).TI=info.InversionTime;
end

%determine sigma of the noise
imshow(fitdata(2).image,[]);
M=roipoly;
close;

for i=1:nfol-2
    D=fitdata(i).image.*M;
    sdyd(i)=std(nonzeros(D));
    maxsd=max(sdyd);
end

clear D;

%conversion factor from standard to rayleigh distribution
CRLBsigma=1.527*maxsd
%CRLBsigma = 150;

%read in data to fitdata struct
for i=1:nfol-2
    if i==1
       imshow(fitdata(2).image,[]);
%when ROI is to be drawn here:
       M=roipoly;
       close;
%when ROI is imported
       %M=logical(maskimage./maskimage);
       %M=ones(256,256);
       %M(130:145,150:180)=1;
       %M=maskimage;

    end
    fitdata(i).image=fitdata(i).image.*M;

end

for k=1:nfol-2
    TI(k)=double(fitdata(k).TI);
end

%initialize&allocate
TI=TI';
yd=zeros(size(TI,1),size(TI,2));

[maxTI,imaxTI]=max(TI);

if ~roibased
    [row,column]=find(fitdata(2).image);
    nrvoxels=size(row,1);
end;

dims=size(fitdata(1).image);
T1=zeros(dims(1),dims(2));
R1=T1;
A=T1;
B=T1;
STDCRT1=T1;
STDCRR1=T1;
STDCRB=T1;
STDCRA=T1;

% CRT1=T1;
% CRR1=T1;
% CRA=T1;
% CRB=T1;
% T1thresh=T1;
% EigVec=zeros(3,3,nrvoxels);
% EigVal=EigVec;

opt = optimset('fminunc');
opt = optimset(opt,'Diagnostics','off','LargeScale','off','gradObj','on','Display','off','MaxIter',80,'Hessian','off','TolFun',1e-19,'Tolx',1e-7);
%s = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,100],'Startpoint',[11500, 11500, 0],'Maxiter',4);
%f = fittype('a-b*exp(-x*c)','options',s);

   if roibased
         yd=double(fitdata(1).image);
        
         for k=2:nfol-2
                yd=[yd double(fitdata(k).image)];
         end
         
         yd=reshape(yd,size(yd,1)*size(yd,2),1);
         TI=repmat(TI,1,size(yd,1)/(nfol-2));
         TI=TI';
         TI=reshape(TI,size(yd,1),1);
         [nonzerox,nonzeroy]=find(yd);
         yd=yd(nonzerox,1);
         TI=TI(nonzerox,1);
         nrvoxels=1;
    end

for i=1:nrvoxels
    disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels(1))]);
        if ~roibased
            if fitdata(1).image(row(i),column(i)) > 0;
                for k=1:nfol-2
                    yd(k)=double(fitdata(k).image(row(i),column(i)));
                end

                [LS,ML] = T1IRComputePar(yd,TI,0,CRLBsigma,opt);
                [CR] = T1IRCramerRao(ML,TI,CRLBsigma);
                
                %[EigVec(:,:,i),EigVal(:,:,i)]=eig(CR);

                A(row(i),column(i))=ML(1,1);
                B(row(i),column(i)) = ML(2,1);
                T1(row(i),column(i)) = 1/ML(3,1);
                R1(row(i),column(i)) = ML(3,1);

                %CRA(row(i),column(i)) = CR(1,1);
                %CRB(row(i),column(i)) = CR(2,2);
                %CRR1(row(i),column(i)) = CR(3,3);
                %CRT1(row(i),column(i)) = CR(3,3)/(ML(3,1)^4);

                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                else
                    STDCRA(row(i),column(i)) = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB(row(i),column(i)) = sqrt(CR(2,2));
                else
                    STDCRB(row(i),column(i)) = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1(row(i),column(i)) = sqrt(CR(3,3));
                    STDCRT1(row(i),column(i)) = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1(row(i),column(i)) = 5;
                    STDCRT1(row(i),column(i)) = 10000;
                end
            else
                A(row(i),column(i))= 0;
                B(row(i),column(i)) = 0;
                T1(row(i),column(i)) = 0;
                R1(row(i),column(i)) = 0;

                %CRA(row(i),column(i)) =0;
                %CRB(row(i),column(i)) = 0;
                %CRR1(row(i),column(i)) = 0;
                %CRT1(row(i),column(i)) = 0;
                STDCRR1(row(i),column(i)) = 0;
                STDCRT1(row(i),column(i)) = 0;
            end  
        else
                [LS,ML]=T1IRComputePar(yd,TI,0,CRLBsigma,opt);
                [CR]=T1IRCramerRao(ML,TI,CRLBsigma);
                A = ML(1,1);
                B = ML(2,1);
                T1 = 1/ML(3,1);
                R1 = ML(3,1);
                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA = sqrt(CR(1,1));
                else
                    STDCRA = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB = sqrt(CR(2,2));
                else
                    STDCRB = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1 = sqrt(CR(3,3));
                    STDCRT1 = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1 = 5;
                    STDCRT1 = 10000;
                end
        end
            
end

%write result in dcm files. R1 values*10000 to get in a sensible range

%lx=nonzeros(EigVec(1,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%ux=nonzeros(EigVec(1,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));
%ly=nonzeros(EigVec(2,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%uy=nonzeros(EigVec(2,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));

if ~roibased
    %mkdir([maskdir,'\',maskfile,'_mappings']);
    dicomwrite(int16(10000*R1),[outputdir,filesep,file,'_R1map','.dcm'],fitdata(1).info);
    dicomwrite(int16(T1),[outputdir,filesep,file,'_T1map','.dcm'],fitdata(1).info);
    %dicomwrite(int16(A),[maskdir,'\',maskfile,'Amap','.dcm'],fitdata(1).info);
    %dicomwrite(int16(B),[maskdir,'\',maskfile,'Bmap','.dcm'],fitdata(1).info);
    dicomwrite(int16(10000*STDCRR1),[outputdir,filesep,file,'_R1errormap','.dcm'],fitdata(1).info);
    dicomwrite(int16(STDCRT1),[outputdir,filesep,file,'_T1errormap','.dcm'],fitdata(1).info);
    imshow(T1,[]);

end

if roibased
    scatter(TI,yd,'.b')
    hold on;
    ti=floor(min(TI)):1:round(max(TI));
    yy=ML(1,1)-ML(2,1)*exp(-ti*ML(3,1));
    plot(ti,yy,'k');
    disp(['Results: value +- standard deviation:'])
    disp(['R1 = ' num2str(ML(3,1)) ' +- ' num2str(sqrt(CR(3,3)))])
    disp(['T1 = ' num2str(1/ML(3,1)) ' +- ' num2str(sqrt(CR(3,3))/(ML(3,1)^2))])
end;

if niftioutput

    T = [0 1 0 0;-1 0 0 0;0 0 -1 0;0 0 0 1]; %useplan(k,2).Tf;
    voxsp = sqrt(sum(T(:,1:end-1).^2,1));

    R1nii = make_nii(R1, voxsp);
    R1nii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    R1nii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    R1nii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    R1nii.hdr.hist.sform_code = 1;
    save_nii(R1nii, [studydir,'\',fitdata(1).folder,'\','R1mapnii'])

    T1nii = make_nii(T1, voxsp);
    T1nii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    T1nii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    T1nii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    T1nii.hdr.hist.sform_code = 1;
    save_nii(T1nii, [studydir,'\',fitdata(1).folder,'\','T1mapnii'])

    Anii = make_nii(A, voxsp);
    Anii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    Anii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    Anii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    Anii.hdr.hist.sform_code = 1;
    save_nii(Anii, [studydir,'\',fitdata(1).folder,'\','Amapnii'])

    Bnii = make_nii(B, voxsp);
    Bnii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    Bnii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    Bnii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    Bnii.hdr.hist.sform_code = 1;
    save_nii(Bnii, [studydir,'\',fitdata(1).folder,'\','Bmapnii'])

    R1errnii = make_nii(STDCRR1, voxsp);
    R1errnii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    R1errnii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    R1errnii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    R1errnii.hdr.hist.sform_code = 1;
    save_nii(R1errnii, [studydir,'\',fitdata(1).folder,'\','R1errormapnii'])

    T1errnii = make_nii(STDCRT1, voxsp);
    T1errnii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
    T1errnii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
    T1errnii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
    T1errnii.hdr.hist.sform_code = 1;
    save_nii(T1errnii, [studydir,'\',fitdata(1).folder,'\','T1errormapnii'])

end

%imshow(T1,[])

%To show an A vs. R2 map with errorbars(either parallel with x and y axis
%or in the direction of the eigenvectors of the covariancematrix) enable one of the following lines: 

%eigenvector errorbaroblique does not work yet for 3 parameter fitting.
%errorbarxyoblique(nonzeros(A)',nonzeros(T1)',lx',ly',ux',uy','w') 
%errorbarxy(nonzeros(B)',nonzeros(T1)',nonzeros(STDCRB)',nonzeros(STDCRT1)',[],[],'w');
%hold on
%scatter(nonzeros(B),nonzeros(T1),'.k')

%disp('value +- standard deviation:')
%disp(['Region: ' num2str(maskfile)])
%disp(['T1 = ' num2str(mean(nonzeros(T1))) ' +- ' num2str(std(nonzeros(T1)))])
%disp(['On ' num2str(nrvoxels(1)-counter) ' voxels'])
%disp(['Mean CR = ' num2str(mean(nonzeros(STDCRT1))) ' +- ' num2str(std(nonzeros(STDCRT1)))])
    
